home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Transactor
/
Transactor_17_1987_Transactor_Publishing.d64
/
n-body simulator
(
.txt
)
< prev
next >
Wrap
Commodore BASIC
|
2023-02-26
|
8KB
|
319 lines
100 rem n-body simulator
110 rem version 6.09
120 rem by richard lucas
130 :
140 rem initialize
150 poke2038,peek(55):poke2039,peek(56)
160 poke 56,62: poke 55,0: clr
170 dim x(40),y(40),z(40),u(40),v(40),w(40),x1(40),y1(40),z1(40)
180 dim m(40),gm(40),e2(7)
190 dim x0(40),y0(40),z0(40),u0(40),v0(40),w0(40),a0(40),b0(40),c0(40)
200 dim ex(7),un(3),un$(3),tu$(3)
210 fori=1to3:readtu$(i),un$(i):next
220 data "days","10^7kilometers-tons-days","days","[193][213]-[197]arth mass-days","seconds"
230 data 1000km-kg-sec
240 gosub 2820
250 g=6.67e-11:sy=1:cf=40:c4=504:c7=7:c8=248:vi=53248:hi=vi+16:c5=255
260 hr=16*1024:o$="epncdls"+chr$(133)+chr$(136)+"xyzuvwmtr[196]":sp=0:dv=8
270 nb=0:dt=10:d2=dt*dt/2:d3=d2*dt/3:d4=d3*dt/4:cb=1
280 fori=0to7:ex(i)=2^(7-i):next
290 fori=0to7:e2(i)=2^i:next
300 un(1)=1000*86400^2/1e10^3
310 un(2)=5.9742e24*86400^2/1.4959789e11^3
320 un(3)=1/1e6^3
330 poke 53280,0:poke 53281,0:print chr$(14)+chr$(8)+chr$(151);
340 b$="":l$="":fori=1to38:b$=b$+chr$(32):l$=l$+chr$(192):nexti
350 b$=chr$(29)+b$+chr$(29)
360 :
370 rem main loop
380 gosub 2460
390 gosub 2620
400 print ">";
410 get a$:ifa$=""then 410
420 fori=1tolen(o$):ifa$=mid$(o$,i,1)then printa$:goto 450
430 next
440 goto 410
450 on i goto 490,510,1970,2900,3070,2110,2250,3190,3230
460 if i=18 then input"[206]ew name of body";n$(cb):goto 370
462 if i=19 then input"[196]rive #";dv:goto 370
470 input"new value";nv
480 on i-9 goto 2390,2400,2410,2420,2430,2440,2380,2450
490 poke55,peek(2038):poke56,peek(2039):clr:end
500 :
510 rem plot trajectories
520 if nb=0 then print "[206]o bodies in current system.":goto 370
530 if sp then 610:rem skip hires for sprites
540 rem set up hires screen
550 poke 56578,peek(56578)or3:rem switch to vic bank 1 (16k-32k)
560 poke 56576,(peek(56576)and252)or2
570 poke 53272,(peek(53272)and15)or128:rem char screen is in 9th k
580 poke 53265,peek(53265)or32:rem turn on hires screen
590 poke 820,0:poke 821,64:poke 822,0:poke 823,96:poke 251,0:sys 49152
600 poke 820,0:poke 821,96:poke 822,231:poke 823,99:poke 251,16:sys 49152
610 ifspthenfori=vitohi:pokei,.:next
620 ifspthen a1=0:fori=1to8:a1=a1 or(-(i<=nb)*2^(i-1)):next:poke vi+21,a1
630 ifspthenpoke 53281,0:printchr$(147)
640 t=0
650 :
660 rem move start parameters to working arrays
670 fori=1tonb
680 x(i)=x0(i)
690 y(i)=y0(i)
700 z(i)=z0(i)
710 u(i)=u0(i)
720 v(i)=v0(i)
730 w(i)=w0(i)
740 next
750 rem compute accel at time dt before start
760 fori=1tonb
770 a0(i)=.:b0(i)=.:c0(i)=.
780 next
790 fori=1tonb
800 ax=.:ay=.:az=.
810 forj=1tonb
820 ifi=jthen910
830 dx=x(j)-x(i)
840 dy=y(j)-y(i)
850 dz=z(j)-z(i)
860 r=sqr(dx*dx+dy*dy+dz*dz)
870 r3=r*r*r/gm(j)
880 ax=ax+dx/r3
890 ay=ay+dy/r3
900 az=az+dz/r3
910 next
920 x1(i)=x(i)-u(i)*dt+ax*d2
930 y1(i)=y(i)-v(i)*dt+ay*d2
940 z1(i)=z(i)-w(i)*dt+az*d2
950 next
960 fori=1tonb
970 ax=.:ay=.:az=.
980 forj=1tonb
990 ifi=jthen1080
1000 dx=x1(j)-x1(i)
1010 dy=y1(j)-y1(i)
1020 dz=z1(j)-z1(i)
1030 r=sqr(dx*dx+dy*dy+dz*dz)
1040 r3=r*r*r/gm(j)
1050 ax=ax+dx/r3
1060 ay=ay+dy/r3
1070 az=az+dz/r3
1080 next
1090 a0(i)=ax
1100 b0(i)=ay
1110 c0(i)=az
1120 next
1130 :
1140 rem calculate new system state
1150 fori=1tonb
1160 a1=.:b1=.:c1=.:a2=.:b2=.:c2=.
1170 forj=1tonb
1180 ifi=jthen1270
1190 dx=x(j)-x(i)
1200 dy=y(j)-y(i)
1210 dz=z(j)-z(i)
1220 r=sqr(dx*dx+dy*dy+dz*dz)
1230 r3=r*r*r/gm(j)
1240 a1=a1+dx/r3
1250 b1=b1+dy/r3
1260 c1=c1+dz/r3
1270 next
1280 j0=(a1-a0(i))/dt
1290 k0=(b1-b0(i))/dt
1300 l0=(c1-c0(i))/dt
1310 x2=x(i)+u(i)*dt+a1*d2+j0*d3
1320 y2=y(i)+v(i)*dt+b1*d2+k0*d3
1330 z2=z(i)+w(i)*dt+c1*d2+l0*d3
1340 forj=1tonb
1350 ifi=jthen1440
1360 dx=x(j)-x2
1370 dy=y(j)-y2
1380 dz=z(j)-z2
1390 r=sqr(dx*dx+dy*dy+dz*dz)
1400 r3=r*r*r/gm(j)
1410 a2=a2+dx/r3
1420 b2=b2+dy/r3
1430 c2=c2+dz/r3
1440 next
1450 j1=(a2-a1)/dt
1460 k1=(b2-b1)/dt
1470 l1=(c2-c1)/dt
1480 m1=(a2-2*a1+a0(i))/(dt*dt)
1490 n1=(b2-2*b1+b0(i))/(dt*dt)
1500 o1=(c2-2*c1+c0(i))/(dt*dt)
1510 x1(i)=x(i)+u(i)*dt+a1*d2+j1*d3+m1*d4
1520 y1(i)=y(i)+v(i)*dt+b1*d2+k1*d3+n1*d4
1530 z1(i)=z(i)+w(i)*dt+c1*d2+l1*d3+o1*d4
1540 u1(i)=u(i)+a1*dt+j1*d2+m1*d3
1550 v1(i)=v(i)+b1*dt+k1*d2+n1*d3
1560 w1(i)=w(i)+c1*dt+l1*d2+o1*d3
1570 a0(i)=a1
1580 b0(i)=b1
1590 c0(i)=c1
1600 next
1610 :
1620 fori=1tonb
1630 x(i)=x1(i)
1640 y(i)=y1(i)
1650 z(i)=z1(i)
1660 u(i)=u1(i)
1670 v(i)=v1(i)
1680 w(i)=w1(i)
1690 x=x(i):ifx(i)<.orx(i)>319then1720
1700 y=y(i):ify(i)<.ory(i)>199then1720
1710 gosub2710
1720 next
1730 t=t+dt
1740 if sp then printchr$(19);t
1750 geta$:ifa$=""then1150
1760 :
1770 rem restore character screen
1780 poke53265,peek(53265)and223
1790 poke 56578,peek(56578)or3
1800 poke 56576,(peek(56576)and252)or3
1810 poke 53272,(peek(53272)and15)or16
1820 poke vi+21,0:rem turn off sprites
1830 poke 53281,0
1840 if a$=chr$(133) then 1860
1850 goto 370
1860 print chr$(17);"[211]toring present system in memory.":gosub 3160
1870 for i=1 to nb
1880 x0(i)=x(i)
1890 y0(i)=y(i)
1900 z0(i)=z(i)
1910 u0(i)=u(i)
1920 v0(i)=v(i)
1930 w0(i)=w(i)
1940 next i
1950 goto 370
1960 :
1970 rem get new system from user
1980 input"[206]umber of bodies";nb
1990 ifnb<1ornb>50then370
2000 for i=1 to nb
2010 cb=i:gosub 2460
2020 print"[206]ame of body"i;:inputn$(i)
2030 iflen(n$(i))>25then2020
2040 print"[205]ass of body"i;
2050 input m(i):gm(i)=g*m(i)*un(sy)
2060 input"[201]nput location in x,y,z form";x0(i),y0(i),z0(i)
2070 input"[201]nput velocity in u,v,w form";u0(i),v0(i),w0(i)
2080 next i
2090 goto 370
2100 :
2110 rem load system description from disk
2120 print "[204]oad system data from disk."
2130 input"[212]ype name of data file";a$
2140 if len(a$)>13 then print"[212]oo long.":goto 2130
2150 open 15,dv,15:open 2,dv,2,"0:"+a$+".nb,s,r"
2160 gosub 3170
2170 if er<>0 then print er$(1);er$(2);er$(3);er$(4):gosub 3160:goto 2230
2180 input#2,sy,nb
2190 fori=1tonb
2200 input#2,n$(i),x0(i),y0(i),z0(i),u0(i),v0(i),w0(i),m(i)
2210 gm(i)=g*m(i)*un(sy)
2220 next
2230 close 2:close 15:cb=1
2240 goto 370
2250 rem save current system to disk
2260 print "[211]ave current system."
2270 input"[212]ype name of file";a$
2280 if len(a$)>13 then print"[206]ame too long.":goto 2250
2290 open 15,dv,15:c$=chr$(13)
2300 open 2,dv,2,"0:"+a$+".nb,s,w"
2310 gosub 3170
2320 if er<>0 then print er$(1);er$(2);er$(3);er$(4):gosub 3160:goto 2230
2330 print#2,sy;c$;nb
2340 fori=1tonb
2350 print#2,n$(i);c$;x0(i);c$;y0(i);c$;z0(i);c$;u0(i);c$;v0(i);c$;w0(i);c$;m(i)
2360 next
2370 close 2:close 15:goto 370
2380 m(cb)=nv:gm(cb)=g*m(cb)*un(sy):goto 370
2390 x0(cb)=nv:goto 370
2400 y0(cb)=nv:goto 370
2410 z0(cb)=nv:goto 370
2420 u0(cb)=nv:goto 370
2430 v0(cb)=nv:goto 370
2440 w0(cb)=nv:goto 370
2450 dt=nv:d2=dt*dt/2:d3=d2*dt/3:d4=d3*dt/4:goto 370
2460 rem display current system values
2470 printchr$(147)" [206]-[194][207][196][217] [211][201][205][213][204][193][212][207][210]"
2480 printchr$(176);l$;chr$(174);
2490 l1$=chr$(221)
2500 printl1$" name: "n$(cb);tab(79);l1$;
2510 printl1$" body #"cb;tab(17)"mass:"m(cb);tab(39)l1$;
2520 printl1$" x:"x0(cb);tab(60)"u:"u0(cb);tab(79);l1$;
2530 printl1$" y:"y0(cb);tab(20)"v:"v0(cb);tab(39);l1$;
2540 printl1$" z:"z0(cb);tab(60)"w:"w0(cb);tab(79);l1$;
2550 printchr$(173);l$;chr$(189)
2560 print "number of bodies:"nb
2570 print "time interval:"dt;tu$(sy)
2580 print "unit system: "un$(sy)
2590 if sp then print "sprite mode";chr$(17):return
2600 print "hires point mode";chr$(17):return
2610 :
2620 rem display menu
2630 print "e[146]xit p[146]lot n[146]ew system"
2640 print "sc[146]ale d[146]isplay l[146]oad [196][146]rive#"
2650 print "s[146]ave f1[146] prev body f7[146] next body"
2660 print "x[146] position y[146] position z[146] position"
2670 print "u[146] velocity v[146] velocity w[146] velocity"
2680 print "m[146]ass t[146]ime r[146]ename"
2690 return
2700 rem plot point on hires screen
2710 if sp=1 then 2750
2720 ml=hr+(yandc8)*cf+(yandc7)+(xandc4)
2730 poke ml,peek(ml)orex(xandc7)
2740 return
2750 if i>8 then return
2760 x=x+24:y=y+50
2770 poke vi+(i-1)*2,xandc5
2780 poke vi+i*2-1,y
2790 ifx>c5thenpokehi,peek(hi)ore2(i-1)
2800 ifx<256thenpokehi,peek(hi)and(c5-e2(i-1))
2810 return
2820 rem ml code for high speed erase
2830 i=49152
2840 readmc:ifmc=256thenreturn
2850 pokei,mc:i=i+1:goto2840
2860 data173,52,3, 133,2, 173,53,3, 133,3
2870 data165,251, 160,0, 166,3
2880 data145,2, 236,55,3, 208,7, 166,2, 236,54,3, 240,9
2890 data230,2, 208,236, 230,3, 76,14,192, 96, 256
2900 print chr$(147);"[211]elect a system of units."
2910 print chr$(17);"1. 1 pixel = 10^7 kilometers"
2920 print " 1 mass = 1000 kilograms"
2930 print " 1 time = 1 day"
2940 print chr$(17);"2. 1 pixel = 1 [193][213] (earth radius)"
2950 print " 1 mass = 1 earth mass"
2960 print " 1 time = 1 day"
2970 print chr$(17);"3. 1 pixel = 1000 kilometers"
2980 print " 1 mass = 1 kilogram"
2990 print " 1 time = 1 second"
3000 print: input"[215]hich system";sy
3010 if sy<1orsy>3then370
3020 fori=1tonb
3030 gm(i)=g*m(i)*un(sy)
3040 next
3050 goto 370
3060 :
3070 rem switch plot systems
3080 if sp=1 then sp=0:goto 370
3090 sp=1
3100 fori=15872to15872+8*64:pokei,.:next:rem blank out sprite images
3110 fori=0to7:poke15872+i*64,224:poke15875+i*64,224:next:rem form dot shape
3120 fori=0to7:poke2040+i,248+i:next:rem set sprite pointers
3130 fori=0to7:pokevi+39+i,i+1:next:rem set sprite colors
3140 poke vi+29,0 :poke vi+23,0 :rem compress sprites
3150 goto 370
3160 forde=1to1500:next:return
3170 input#15,er$(1),er$(2),er$(3),er$(4)
3180 er=val(er$(1)):return
3190 cb=cb-1:if cb<1 then cb=nb
3200 printchr$(19)chr$(17)chr$(17);b$;b$;b$;b$;b$:printchr$(19)
3210 gosub 2480
3220 fori=1to6:printchr$(17);:nexti:printchr$(29);:goto410
3230 cb=cb+1:if cb>nb then cb=1
3240 goto 3200
3250 poke53272,peek(53272)and247
3260 poke53265,peek(53265)and223